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ABSTRACT 


This two-part report is concerned with the development of a general framework 
for the implicit time-stepping integrators for the flow and evolution equations in 
generalized viscoplastic models. The primary goal is to present a complete theoretical 
formulation, and to address in detail the algorithmic and numerical analysis aspects 
involved in its finite element implementation, as well as to critically assess the numerical 
performance of the developed schemes in a comprehensive set of test cases. On the 
theoretical side, the general framework is developed on the basis of the unconditionally- 
stable, backward-Euler difference scheme as a starting point. Its mathematical structure is 
of sufficient generality to allow a unified treatment of different classes of viscoplastic 
models with internal variables. In particular, two specific models of this type, which are 
representatives of the present state-of-art in metal viscoplasticity, are considered in 
applications reported here; i.e., fully associative (GVIPS) and non-associative (NAV) 
models. The matrix forms developed for both these models are directly applicable for both 
initiall y isotropic and anisotropic materials, in general (three-dimensional) situations as 
well as subspace applications (i.e., plane stress/strain, axisymmetric, generalized plane 
stress in shells). On the computational side, issues related to efficiency and robustness are 
emphasized in developing the (local) iterative algorithm. In particular, closed-form 
expressions for residual vectors and (consistent) material tangent stiffiiess arrays are given 
explicitly for both GVIPS and NAV models, with their maximum sizes “optimized” to 
depend only on the number of independent stress components (but independent of the 
number of viscoplastic internal state parameters). Significant robustness of the local 
iterative solution is provided by complementing the basic Newton-Raphson scheme with a 
line-search strategy for convergence. In the present second part of the report, we focus on 
the specific details of the numerical schemes, and associated computer algorithms, for the 
finite-element implementation of GVIPS and NAV models. 
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Robust Integration Schemes for Generalized Viscoplasticity 
with Internal-State Variables; Part II Algorithmic 
Developments and Implementation 

1. Introduction 

The scope of the work in this report focuses on the implementation and algorithmic 
developments of two classes of viscoplastic models: GVIPS (folly-associative) and NAV 
(nonassociative) based on the theory discussed in Part I [41] of the report. In the computer 
implementation of a viscoplastic model, the computational algorithm is the key ingredient. Over 
the past years considerable research effort has been devoted to the development of 
computational algorithms [1-9, 12, 16-50], As discussed in Part I [41] of the report, initially, 
simple explicit integration schemes were predominate in finite element applications because of 
their ease in implementation, and because they do not require evaluating and inverting a 
Jacobian matrix. However, explicit integrators may not be efficient. That is, too many iteration 
steps may be required and convergence stability can not be guaranteed [32, 36, 50]. As a 
result, several alternative approaches have been used, for example, Gear’s multi-step method 
[14] and Walker’s asymptotic method [47], Note that every integration scheme has its own 
particular application domain and is very problem dependent. For practical problems, it is 
desirable to use an unconditionally stable integration scheme in order to obtain an accurate 
solution. Recent work has clearly emphasized the use of implicit integration methods [4-6, 12, 
16, 19-20, 23-24, 29, 31, 35, 37-38, 40, 43-44, 46, 48-49] in view of their superior stability 
and convergence properties [18, 34]. From the standpoint of practical applications, the one- 
step, fully implicit, backward Euler scheme is presently one of the most widely used integrators 
[6, 16, 18-20, 23, 31, 34, 35, 38, 39, 43, 44, 48, 49], For the computationally intensive 
viscoplastic applications, found in typical finite element analysis, implicit backward Euler 
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integration methods have become the proven standard for the numerical integration of the 
viscoplastic rate equations [IS, 40]. Details on the implementation of both classes of 
viscoplastic models (NAV, GVIPS) will be discussed in this report. 

Based on the fully implicit , backward Euler scheme, the corresponding algorithmic 
(consistent) tangent stillness arrays are derived from the integration rule, which are important 
for finite element solutions using (global) Newton-Raphson iterative methods. Explicit, dosed- 
form, expressions for state variables and consistent tangent stiffness matrixes are given for the 
general three-dimensional (3D) situations, as well as for their direct modifications for the 
effident treatment of two-dimensional (2D) or subspace applications; i.e., "generalized" plane 
stress states in shells. For these latter 2D problems, this simply amounts to "appropriate" 
reduction in the dimensions of the arrays involved, which is known to be more effective than 
alternative implementations [e g., 12, 49], in which indirect (iterative) methods are used to 
handle the zero- stress constraints. 

Beginning with the second section, and for the remainder of the report, all equations and 
expressions are shown using concise matrix notations for the integrated stress and internal 
stress fields. The contracted (Voigt) representations in vector forms for the components of the 
corresponding symmetric, second-order, tensor are utilized and the appropriate dimensions for 
the space of these vectors are defined explidtly; i.e., a six-dimensional space for 3D continuum 
problems, etc. With a slight abuse in notations to indicate the vector-matrix representations, 
matrices are defined by under-curved symbols, and vectors by underlined symbols, and ® is 
used to indicate the tensor product of two vectors. Considering vectors a, b and c of 
dimension (6x1) in the 3D case, we define the matrix (a 0 b)c=(b • c)a, where signifies the 
inner (scalar) product of two vectors, and ": M for a double contraction C:D (or Qju Du) and 
for tensor product u 0 w (or Ui Wj). Note that ef_ indicates the inelastic strain. 

Section 2 deals with the implementation of the GVIPS class of models, as described in Part 
I [41], and similarly section 3 explains the implementation of the NAV class models. Although 
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the implicit backward Euler integration scheme is stable, for highly-nonlinear viscoplastic 
models, convergence may be difficult for certain problems. As a result, line search, a numerical 
technique based cm optimization theory, is utilized to guarantee convergence and improve 
effi ciency of the integrator. Details of the line search method are discussed in section 4. 


1 . 1 General Form of Newton Iterative Scheme 


The backward-Euler scheme is based on the equation 


~ 5ii + A2* +1 


( 1 . 1 ) 


where Zn is state variables, and is the increment of the state variables, n is the step 
counter. For the GVIPS models. 


2 - =1 


a 


n) 


( 1 . 2 ) 


for the NAV models, 


2 .= 




a. 


<Xj_ 


(1.3) 


the expression for the state variables’ increment is 


AS„, = (a) 


(1.4) 
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where K z is iterative Jacobi matrix of state variables, R n+l is the residual function of state 
variables. The specific forms of K z and R„ +l for GVIPS and NAY are given subsequently. 


2. Implicit Integration Scheme: GVIPS Class 

The details of the integration algorithm, using the fully-implicit Euler method, for the 
specific GVIPS model (described in Part I [41]) are given below For convenience in 
implementation, all equations of this GVIPS model are written in a vector-matrix format This 
model serves as an example for discussing the implicit integration scheme. 

2.1. The General Form 

The hyperelastic response of the model is assumed to be linear, i.e., the Cauchy (true) 
stress components a are given by 


& = C‘(c- £ f ) ; e = e‘+e p 

(2.1) 

The governing equations for GVIPS are: 


e p = f{F)T ; r=M(a-a) 

(2.2) 

a = Zr'f*^ j ; n = Ma 

(2.3) 


where, 

F = -(a-a)M(a-a)/Kf-l (2.4a) 

G = -a MalK) 

2 ~ ~ ~ ' 


(2.4b) 



5 


in which the symmetric matrix M is a function of the flier direction and is rewritten as 

M = P-ZQ-\£R (2-5) 

~ Zr 

where P, Q, and R are symmetric matrices as defined in Part I [41]. 

In eq. (2.3), 

L" 1 = h 

In the above equations, £ h, and y are material functions which are defined as follows. 


f{F) = ri(2nK t 4F7i) 

(2.6b) 

h = H/G p ; y = RG m-p / (k,Vg) 

(2.6c) 


with, n, m, ti, Kt, H, P, and R denoting material dependent constants. In addition, the 
“generalized inverse” of M is defined as follows: 

Z m = M' 1 (2.7a) 

withM directly evaluated from eq. (2.5) for the case of plane-stress continuum (with dimension 

3x3) and generalized plane stress in shells (with dimension 5 x 5), but for the three- 
dimensional case, P (6x6 dimension) in eq. (2.5) needs to be replaced by 


Z m + 


h' 


h(l + 2p) 


a®a 


(2.6a) 


P = diag[ 1,1,1, 2,2,2] 


(2.7b) 
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with its appropriate direct reduction to (4 x 4) for the axisymmetric/plane-strain continuum 
problems, that is, 

P = diag[ UU] (2.7c) 

In the above, diag [•] indicates a diagonal matrix with entries [-]. 

The basic problem considered here is as follows. Consider a typical time step 
where the state variables at t„ are known; i.e., {a*, ctn, e£}. Let A t=t n *i~ t„, and 

Ae = Ate =e n+1 -e n be the given increment in the total strain field. Based upon the quantities 
at it is required to update to time i„+j the above fields {cvi, <vi, e^-i} in a manna- 

consistent with the governing equations; i.e., eqs. (2.2, 2.3) for £ 'and a, and the rate form of 

eq. (2.1) for d, i.e., 

& = C(e-e p ) ; e = e'+e p (2.8a) 

the matrix form of eq. (2.3) may be written as, 

s=f j/r- 

Using the implicit Euler scheme, these update formulas are given as 

CLn+i =&,+£•[ Ae- AtfT^ ] (2.9a) 

On+x = Z„+At L~' [/£*, - j (2.9b) 

where 

f=f(a*+i, ctn+i) ; h = h( aLn+ i) ; Y = Y (a^i) 


Y ~ 

— TV 

h~\ 


(2.8b) 


( 2 . 10 ) 
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2.2 Iterative Solution for Updated Fields 


To determine the updated values for the state variables (3*1, 05^1), a local iterative 
solution (i.e., distinct from global equilibrium iterations in finite elements) is needed based on 
the system of equations in Eqs. (2.1 and 2.3). To this end, we apply the Newton-Raphson 
scheme by first forming the residual vectors 

B, = 2». - [e. + £“(Af - A0T)] 

(2.11a) 


(2.11b) 

A truncated Taylor's expansion for the R vectors about the last updated state that yields, for a 

typical iteration k—>k+\ : 


k+\ k . a k 

(2.12a) 

®«+l = — n+1 + A — n+1 

(2.12b) 

where 



(2.13a) 


(2.13b) 

in which J is the iteration Jacobian matrix, and the following definitions are 

introduced: 

J = P, - P, P, P, = C*"' + P, P, j 

(2.14a) 

P, = Z m ~'+P 2 + AtyM+f"Aty'--^-(Aty 

(2.14b) 
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P 4 =At(fM+f'r®£) ; P 2 =hP 4 

(2.14c) 

P 3 = C'~ l + At{fM+ /T®zz) 

(2. 14d) 

K = (m- - a „ ) - At(hfT - yn) 

(2. 14e) 


and the primes indicate "scaled" derivatives as given below (see Eqs. 3.30 in Part I [41]): 


J k 2 ,dF 


7 dh 1 df 

JldG '' T ~~kjdG 


(2.15) 


Note that all the arrays and functions, e.g. L, F, Jt, fj h etc. in Eqs. (2.13) to (2.14) are 

evaluated based on the last updated state (o^i, a^i). Note also that the state (a* a*) 
corresponds to the converged global solution at the last time step; i.e., they are kept constant 
during the local iteration sequence in Eq. (2. 12). 


2.3. Iterative Solution — Consistent Tangent Stillness 

In keeping with the underlying fully implicit integration scheme, one can 
straightforwardly proceed to differentiate the updated stress field, cvi, to obtain the 
viscoplastic-moduli or material tangent stiffness matrix, C ev , for use in the finite element 
stiffness calculations. The derivations will in fact lead to 

C" = <?ct„ +1 fdAe ; do = C" de (2.16a) 

C = J' 1 = [ C*" 1 + P 4 - hP 4 7> 7>] (2. 16b) 


Clearly, C* and P x , P 2 , P 3 , P 4 (see Eq. 2.14) are symmetric. 
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3. Implicit Integration Scheme: NAV Class 

In this section, Freed’s viscoplastic model [13] is used as an example to discuss the 
implementation scheme of NAV. The discussion here is limited to the case of an isotropic 
material under isothermal conditions. A flowchart of the integration scheme is also 
included in this section. 


3.1 General Form 

The general form of the integration algorithm described in section 2 may also be used 
for a nonassociative model. Here the framework of section 2 is utilized to recast the 
equations of NAV, (given in Part I [41]) and develop the following numerical algorithm. 
Assume linear hyperelastic response of the model, the Cauchy (true) stress components a 
are given by 


a = C‘(£-£ p ) ; §_ = §L + §L (3.1) 

and C* is the elastic material matrix. According to eqs. (3.42, 3.47 and 3.49) in Part I 

[41], the governing equations for NAV are 

eL = f{J,D) £ ; T = M(a-a) (3.2) 

a^ = 2Z(HJT-g t ^j ; r,=A/a, (3.3a) 

a l = 2Z(H l /T-g,7^j ; ^ = A/a^ (33b) 

b = qj -q D (34) 

where the material functions f, g,, gi, qj, qD are defined as follows. 


f = 9A 


2V7 


Sirth* 


g t =^-3ASinh n 
* s 2 L. 


D 

(VT| 


V 


(3.5a) 

(3.5b) 



10 


gl = 


A 

2L, 


SASinh" 



(3.5c) 


D-D, 

m 

/ 

8C 

SASinh” 

V 

»4 Z, - D 0 

\ SC J. 



J = ^(a - a) T M{q_ - a) 


(3.5d) 


(3.5e) 

(3.5f) 


Note that S, A, H* Hi, m, n, D 0 , C and 5 are material constants, and ho is a function of the 
drag stress D (see eq. 3.53 in Part I [41]). Note also that due to the assumption of an 
isotropic material, £, C, found in the expression for M , (eq. 2.5) are set to zero, thus, 


M=P. In this simpler case, following a similar procedure as in section 2, we form the 


matrix Z, termed the generalized inverse of the deviatoric projection P, for different 
spaces; e.g., in three-dimensions, we have 

Z = (p+i<S®£) (3.6a) 

where diag [.] indicates a diagonal matrix with entries in [.]. Also, for generalized plane 
stress (plates/shells), we have: 


Z = P ' ; P = (5x5) (fromP ineq. (3.24a) in Part I [41] ) (3 .6b) 


Again, the objective in solving the incremental problem is to find cr n+1 ,a n+1 and 
D B+1 at time t„ + i based upon the converged values, <r n ,a n ,D n and Ae n at time t„ for the 
given At and Ag whereas n is the step counter. From the rate form of, &, eq. (3 .1), the 
evolution equations (3 .2-3.4) a s , a, and D, and by using the implicit Euler scheme, the 
following expressions are derived, 
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2^1 =£„+C‘ Ae- Atfr„, 

(3.7) 

=£ i .+ 2 A/z(//jr„, 

(3.8a) 

Hi.., = Hi. +24 'z(/f,/r„ l i„,) 

(3 8b) 

A»+i = -At + — Qd) 

(3.9) 

For a local Newton-Raphson scheme to update state variables, it is necessary to form 

residual vectors in terms of the variables, i.e., 


R„ = 2^1 -*.-<?&£+ WZ H+X 

(3,10) 

Ea, =^i„ +1 -g,.-2A tHJ zr^ +2A tg,Z7c i ^ x 

(311a) 

— a i = ^Lin+i — - — zr^, + 2A/£, z^-, 

(3.11b) 

jg = 'Rg, + -Rg, 

(3.11c) 

n 

+ 

1 

i 

Er 

? 

1 

£ 

'w' 

(3-12) 

The local iterative update expressions are based on a truncated Taylor’s 
R vectors about the last updated state, at a typical iteration k — > k+1 

expansion for the 

*+l k A k 

&n + 1 ” Q-n+\ + ^£*+1 

(3.13) 

k + 1 * , * k 

a s =a s +Aa t 

— £*+l — iw+1 — iiH-1 

(3.14a) 

*+l k . A k 

a j = a, +Aa x 

— L/»+i — L#i+i — -»+i 

(3.14b) 

k + 1 ^ *+l , _ k+\ 

a-' =£<„, + £i„, 

(3.14c) 

D‘J = D‘ M +Al& 

(3.15) 


where expressions for the various increments shown above are as follows: 

AgH, = - JT 1 [RAB - PB PC' • RDC) 


(3.16) 
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A a. k ^-PC-^RDC+PD Aal,) 

(3.17a) 

A<,-Vk, -v„ +[-P ll +P 3 ,+ P i , 

(3.17b) 

^A»+i — ^9 ( — S s r ■ Acr^, ^ S r-Aa, ri + s 5 F- 

(3.18) 

In the above, the scalar quantities, S 1 -S 9 , relative to the differentials of the material 
functions in eqs. (3.5), are given as 

- At f D S is — ^AtH x f D Sjj — 2At H x f D 

(3.19a) 

% ^At g sD s 4 = 2 At g lD s 5 = Atqjj s 6 = 0 

(3.19b) 

&j — Atq JD % — Atq DD — 

1 - 57 -* 

(3.19c) 


Where, the second subscript in the primed quantities indicates the variable with respect to 
which the differentiation is performed. The vectors used in eqs (3.16, 3.17 and 3.18) are 
relative to the residual function of state variables and are defined as: 


v* = s,s 9 R d c r 

(3.20a) 

V a, = ‘V^,d(‘ S 2i ~ ZXj) 

(3.20b) 

^ = s 9 R D (s 2l Zr-s 4 Zn L ) ^ = v^ + v^ 

(3.20c) 

+ -v 0 J 

(3.21a) 

RDC=R., -V'. -(p„-P„-P„ V'K -V.) 

(3.21b) 


Among the matrices defined in eqs. (3.16, 3. 17 and 3.18), K is the iterative Jacobi matrix, 
and is defined as 
K=PA-PBPC~ X PD 


(3.22) 
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where, 

PA = UP ,- /;-(V P s ) /),- PB = 1-PA (3.23a) 

^ = Pp = I+P^-PC (3.23b) 

in which matrices P\- P-, result from state variables and their differentials, they are defined 
as 

i> =AtC(fM+ fjr<S>r) (3.24a) 

P 2s = 2 AtH, z(/ M+ fjL®r) (3.24b) 

P 2l = 2 At H, Z( / M+ f'j r®r) (3.24c) 

P 3s = 2 J/Z(^X 0£) P 3 , = 2/fr Z( ft ^0 £) (3.24d) 

P 4s = 2Atg s ZM P 4l =2Atg,ZM (3.24e) 

P 5 = W9 rr®/: (3.240 

P 6 , = s i s 9 (s 2s zr®r-s 3 zx 1 ®r) p 6t =s 5 s 9 (s 2l zr®r-s 4 Z7T l ®r) (3.24 g ) 
P 7 , = I+P 2 ,+P 4 ,-P 3 -P 6 . p v=l +p v +p «- p i‘~ p ‘‘ ( 3 24h ) 


Note that all the arrays and functions in eqs. (3.16-3.24) are evaluated based on the last 
updated state (tr* «/ , «,* and £>* +1 ) which indicate implicit integration scheme. 

Finally, the consistent tangent stiffness may be derived directly from integration shceme 
and is defined as 

(T = K~'C (3.25) 
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Note that K (defined in eq. 3.22) is no longer symmetric, thus leading to unsymmetric C ev , 

which is a major difference between the NAV and GVIPS (whose consistent tangent 
stiffness derived in section 2.3 is symmetric), with significant implication regarding the 
numerical implementation on both the local- and global- stiffness levels. 

3.2 Integration Scheme 

As shown above, GVIPS and NAV have the same implementation format, thus, NAV 
was chosen as the example for presenting the details of the algorithm. The flowchart is 
included in Appendix I. The box labeled Global Iteration represents the calculation of the 
structural stiffness and residual force at the global level. The global structural stiffness is 
based on the current local material stiffnesses that vary with the level of inelasticity 
involved. The purpose of the material model is to integrate the rate form equations (flow 
and evolution laws) of the viscoplastic model over the step size At and obtain the current 
material stiffness required for the update of the structural stiffness and the increment of 
stress A a at the material level. In the flowchart of Appendix I, the subroutines with the 
ml 3 prefix are those necessary to perform the implicit integration of NAV. The main 
subroutine, ml3nrS, is a driver subroutine for the implementation of NAV. It calculates 
increments of stress and internal variables, updates them and checks convergence. Finally, 
it passes the viscoplastic stiffness matrix, and, converged stress to the global level 
calculation. Eqs. (3.13-3.21, 3.25) are involved in this subroutine. Subroutine ml3heav5 
calculates the pretinent material scalar functions (eq. 3.5) and their differentials. 
Subroutine ml3err deals with calculation of the residual functions of stress and internal 
variables during the local iterations. Eqs. (3.10-3.12) are involved in this subroutine. 
Subroutine ml3k calculates iteration Jacobi matrices of stress and internal variables. Eqs. 
(3.22-3.24) are implemented in this subroutine. Subroutine ml31ines performs line search 
algorithm discussed in section 4. Other subroutines used in ml3nr5 are also introduced in 
Appendix I. 
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4 Implementation of Line Search 
4.1 Introduction 

Although the implicit scheme described above is unconditionally stable, its 
successful application still requires proper selection of the size of the steps utilized. In 
this regard two factors are important: (i) accuracy; and (ii) convergence of the local 
iterations. A simple time-subincrementing strategy was found to be effective in 
obtaining accurate results especially when dealing with regions of discontinuity in the 
state space. However, this was found to be insufficient to obtain a computationally 
efficient solution for a highly nonlinear problem, such as viscoplasticity. When a large 
time-step size is chosen, too many subincrements are needed, which leads to 
inefficiency. Thus a more sophisticated solution procedure, namely, a line search 
algorithm, is required to produce an effective, robust solution algorithm. 

The line search technique is an important feature of most numerical techniques 
for unconstrained optimization and can be used with a wide range of iterative solution 
procedures such as full Newton-Raphson iteration. It’s well known that classical 
Newton-Raphson is fast and stable only when the trial solution is close to the 
converged value. For nonlinear problems, the trial solution is usually far away from 
the real solution, thus full step size of iterative increment vector may cause either 
wrong direction or out-of-range updated value. The purpose of the line search 
algorithm is to guide the solution towards convergence, especially when convergence 
becomes more and more difficult, e.g. using excessive number of iterations, 
oscillations in residuals, stresses, and displacement norms, and searches for a scalar 
multiplier that adjusts the amount of the iterative increment vector to be updated 
within each iteration. 
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Line search methods were utilized to solve complex nonlinear problems by 
many researchers. Crisfield [10, 11] applied it to arc-length algorithms to solve 
concrete cracking problems. In elasto-plastic analyses, Simo and Taylor [43] suggest 
that the line searches be incorporated with a consistent tangent stiffness. In [43], Simo 
claims that the use of line search is essential for robust performance of Newton’s 
method. Caddemi and Martin [6] have also demonstrated that in elasto-plastic analysis 
convergence is not guaranteed unless line search is used. 

Among all of the research work mentioned above, the line search method was 
used to optimize the system of the global equations and minimize the out-of-balance 
force. Actually, the concept of line search may be applied at either the global 
(structural) iteration level or at the local (constitutive) iteration level. At the global 
(structural) level, the concept of the line search algorithm pertains to minimizing the 
total potential energy, that is, the work done by the residual force due to the iterative 
displacement. On the local (constitutive) level, it adjusts the suitable increment step of 
stress and internal variables to guarantee the convergence of the stress and internal 
variables at material points, and does not involve global iteration. In this report, the 
line search algorithm was applied to the local level. 

4.2 Line Search Strategy 

The objective of the line search is to optimize the solution system and minimize 
the out-of-balance entropy. Thus, a criterion is needed to judge whether or not a 
given iterative solution is better than a previous one. This criterion takes the form of 
an objective function or cost function. In the following, some basic concepts are 
described that are applicable to unconstrained numerical optimization methods. 

Instead of using classical Newton iteration scheme in eq. (1.1), the following 
iterative procedure was utilized: 


k=0,l,2 


(4.1) 
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In this equation, the superscript k represents the iteration number, Z° is any starting 
value, and A£ k represents an increment value. The iterative scheme described in 
eq.(4.1) is continued until optimality conditions are satisfied, or an acceptable new 
value is obtained. Here, £ k represents the state variable vector associated with a 
viscoplastic constitutive model algorithm. For the GVIPS model, 

( ^ 

E k = (4.2a) 

U J 


and for the NAV model, 

k \ 

a 

£.* (4.2b) 

« i 

D k ) 

Define a scalar function f as the objective function to be optimized as follows: 

= ; df = Rd'L (4.3) 



where, R is the residual function of the state variables. Symbol indicates dot 
product of two vectors. Again for the GVIPS model. 



and for the NAV model, 



(4.4b) 


For unconstrained problems, the residual comes from the difference between the 
current value and the previous converged value. The calculations for AS depend on 
the residual function and its derivatives at the previous step. In order to reach a 
minimum point for the objective function f(S), a suitable value of the scalar r| must be 
found, such that 



18 


f(S k+1 ) = f(I k + iiAS k ) 


(4.5) 


is a minimum. 


Assume AZ k is known at iteration k. According to the stationary condition: 
df(z k+l ) di k+i * 

V /_ V > _ J \- (4.6) 


dr} cFL dr\ ffL 

From eq. (4.3), eq(4.6) becomes: 

R k+1 • AI k = 0 

Now consider, 

s(l) = E(n) ■ A£ k 

in which E k and A£ k are fixed, and E and s are functions of r) When r| = 0 , 
so - s(n=0 ) - A£ k ■ E{ ri=0 ) = A£ k • Ro 


(4.7) 

(4.8) 

(4.9) 


Ro is the residual function of stress and the internal variables at the end of the 
previous iteration. According to optimization theory, the best (optimum) solution is 
s(rj) = 0, but numerically, this is not realistically possible and it is inefficient to try and 
achieve this objective. In practice, a ‘slack’ line search is used, see Fig 1, wherein the 
objective is to make the modulus of s(q) small in comparison to the modulus of So, i.e. 



(4.10) 


where Pi, is the ‘line-search tolerance’. Based on past research experience, a suitable 
value for Pu is on the order of 0.8 [1 1], 


4.3 Search for r| 

In the line search procedure, a key step is the search for an optimum r\ for 
which the requirement ( eq. (4.8) ) is met. Crisfield [11] suggests using a simple linear 
interpolation and extrapolation. Usually no more than two searches for t| are 
necessary. Initially, s 0 /s 0 = 1 at t|o = 0, and the line search is started using rj=l, which 
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corresponds to the basic Newton-Raphson iteration. From eq.(4.8), Si = s(r|=l), thus 
the initial search for a new “improved” value for r\ (denoted as ry in Fig. 2) 
completely depends on the value of Sj/s,, (see point “1” in Fig. 2). Figures (2-3) show 
four possibilities, but the case most frequently encountered is that shown in fig. 2. 

Because of the stiff behavior of the viscoplastic model, A Z may be very large 
especially if a large time-step size is chosen. Consequently, the residual function may 
also be very large, i.e. | Si | may be very large, which leads to a very small rii after 
linear interpolation. Thus, a double interpolation is used by means of ^/sq (denoted as 
point “2” in Fig. 2) and either s 0 /s 0 (Fig. 2a) or s,/so (Fig. 2b) to obtain an updated 
value r )2 which is more accurate than T|i. Note that a minimum of 0.01 is set for ty 

Consider Figure 3 a, which is a possible extrapolation case. When si/so is near 
1.0, extrapolation could result in a large rj which could lead to excessive iterations. 
Even if t| is not very large, as found in the present research a value of r\ which is less 
than 1 is used to speed iteration. For example, Crisfield [11] suggests a maximum 
value of 10.0 for rj when considering concrete cracking analysis. In a similar fashion, 
negative extrapolation of Fig. 3b is not used in the model implementation. 
Theoretically, this could not happen because of the hardening behavior of 
viscoplasticity, only a softening problem may behave like this. Thus, if either of the 
cases described in Figure 3 occur, rj is assumed to be 1 which is the regular Newton- 
Raphson iterative method. 

Numerical tests have shown (refer to Part I [41]) that the use of line search 
leads to significant improvement in the convergence and computational efficiency with 
regards to CPU time and iterations. It has also been shown that line search helps 
convergence greatly for the nonassociative viscoplastic model such that subincrements 
are not required. The only additional computational effort, in comparison with a 
formulation without line searches, would be the calculation of the inner product (eq. 
4.8) which is considered to be almost negligible. 











Ratio 
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Figure 2 Interpolation Schematic 
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Figure 3 Extrapolation Schematic 
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5. Summary and Conclusions 

The general computational framework using a fully implicit backward Euler 
integration method has been presented and has been shown [41] to be successful for both 
the Generalized Viscoplasticity with Potential Structure (GVIPS) and Non-associative 
Viscoplastic models (NAV). Original equations for the NAV model developed by Freed 
are recast into a matrix format similar to that used for GVIPS in order to facilitate the 
model’s implementation into the newly developed framework. Only these Newton iterative 
schemes developed herein will provide, uniformly-valid, convergent, robust integration for 
both GVIPS and NAV. The algorithm was written in a concise matrix format so as to 
provide sufficient generality and is automatically valid for both isotropic and anisotropic 
cases in either full space or subspaces. 

The “slack” method, which is the form of line search used in the present algorithm, 
enables GVIPS and NAV to converge stably, with sufficient accuracy, and significantly 
improved efficiency. Convergence may be achieved for large load steps, whereas 
traditional fixed stepping without line search will fail. Thus, the proposed computational 
framework makes the solution of realistic nonlinear finite element analysis problems 
possible. The only additional computational effort of line search method in comparison 
with a formulation without line search, would be the calculation of the inner product, 
which is believed to be negligible as compared to the total calculations involved in a 
nonlinear finite element analysis. 
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Appendix I: Flowchart for Implicit Integration Algorithm 



ml3nr5 Driver for performing implicit integration of rate-dependent 

equations 

ml3heav5 — Calculate step or smoothing functions 

ml 3 err Calculate error functions during local iteration 

ml 3k Calculate Jacobi matrix during local iteration 

ml31ines Perform line search algorithm 

m 13 j Calculate invariance, i.e., eq. (3.6f) 

ml3maxj Calculate maximum invariance 









Box 1 

Iterative procedure (subroutine m!3nr5 ) 


|l. Given g^o^a^-D^Ae, At. 

k=0, all = ;<*,° B+1 = a, _ ;a,^ = ;£>1, = A, 


L»+l 


|2. Evaluate material function and differentiation. ( subroutine ml3heav5 ) 
£ & qj> Qd> ft> , 

|3. Evaluate residual function ( subroutine ml3err ) 

Ra(5li) , » R o(Dj 

— " +1 = |^2-’ ^a, > A, > R D I 
|4. Evaluate Jacobi matrix K ( subroutine ml3k ) 


|5. Evaluate Ag* +1 , Aa , * +] , Aa , * +] , A D k 
A&ri Ao* + , , Aa,^ , Aa, * , , AD* } T 


L/i+1 


|6. evaluate current value. 

a, M =a, k +Aa / ; .a, k+l = a, k + Aa, k 

— liH-1 — iiH-1 — -w+1 — -n+1 — -n+ 1 — -«■ 


-«+l 




\l. Evaluate residual function based on current value. ( subroutine ml3err ) 

IUO , ^(a,“) , £j(a£;) , J»„(d£S) 

*1*; = 
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Box 1 continued: 





Box 2. 

Line search scheme ( subroutine ml31ines ) 


1. Initialize 

a sL .sIL. = i,tu = o,^ = jSS = rL 

2. Bracketing 

IF sign(A2* +1 • i? £ Jx sign(AI* +1 • <0 then 

interplotating, find rj. 
evaluate = I* +1 + r|AZ* +I 

evaluate residual function based on Zjj +] ( subroutine ml3err ) 

R^ — |^CT > Rg, > Rg, 1 Rp | 

IF sign|A2* +1 • -R T) ) x sign(A2* +1 *R^j <0 then 

interplotaing new rj between rjb and r\ 
else 

interplotaing new r| between r]> and ti 

endif 

ELSE 

il=l 

ENDIF 

IF(r| < Min. ) r|=Min. 

3. Update 

5!L + !=^i+hA2L, 

4. Return. 
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